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Abstract 

We present the first numerical solutions of a kinetic theory description of self- 
sustained current oscillations in n-doped semiconductor superlattices. The governing 
equation is a single-miniband Boltzmann-Poisson transport equation with a BGK 
(Bhatnagar-Gross-Krook) collision term. Appropriate boundary conditions for the 
distribution function describe electron injection in the contact regions. These con- 
ditions seamlessly become Ohm's law at the injecting contact and the zero charge 
boundary condition at the receiving contact when integrated over the wave vector. 
The time-dependent model is numerically solved for the distribution function by us- 
ing the deterministic Weighted Particle Method. Numerical simulations are used to 
ascertain the convergence of the method. The numerical results confirm the validity 
of the Chapman-Enskog perturbation method used previously to derive generalized 
drift-diffusion equations for high electric fields because they agree very well with 
numerical solutions thereof. 
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1 Introduction 



When non-interacting electrons in the conduction band of a material are sub- 
ject to a constant electric field E, their positions should oscillate with a fre- 
quency proportional to the electric field, ojb = eEl/h, where — e < 0, h and 
/ are the charge of the electron, the Planck constant and the crystal period. 
These coherent Bloch oscillations (BO) and the associated current were pre- 
dicted by Zener in 1934 |21| . Scattering limits the observability of BO: to 
observe them, their period should be smaller than the scattering time r, so 
that E > h/ (eh). In standard materials, the fields required to observe BO are 
too large and therefore damped Bloch oscillations were not found until 1992 in 
experiments with undoped semiconductor superlattices which have much 
larger periods than natural crystals. 

Semiconductor superlattices are artificial one-dimensional crystals formed by 
epitaxial growth of layers belonging to two different semiconductors that have 
similar lattice constants [3J. They were synthesized following Esaki and Tsu's 
idea that these artificial crystals would be useful to realize BO or related high 
frequency oscillations [5] . The difference in the energy gaps of the component 
semiconductors makes the conduction band of the superlattice to be a periodic 
succession of barriers and wells with typical periods of several nanometers. The 
electronic spectrum of a superlattice (SL) consists of a succession of minibands 
and minigaps generated by its periodicity. Tayloring the size of barriers and 
wells and the negative doping density of the latter, it is possible to achieve SLs 
with wide minibands and to populate only the lowest one. Electrons moving in 
this miniband have energies that are periodic functions of their wave numbers 
and are scattered by phonons, impurities and other electrons. When an appro- 
priate dc voltage is held between the ends of one such SL with finitely many 
periods, it is possible to obtain high-frequency self-sustained oscillations of the 
current through the structure [JJ. These oscillations are caused by repeated 
formation of electric field pulses at the injecting contact of the SL that move 
forward and disappear at the receiving contact. Thus they are transit-time os- 
cillations whose frequency is inversely proportional to the SL length: they are 
similar to the Gunn effect in bulk semiconductors [TB] and are different from 
BO. These Gunn-type oscillations have been observed in experiments with 
GaAs/AlAs SL (and with other SL based on III-V semiconductors) since 1996 
and are the basis of fast oscillator devices [13]. The connection between the 
existence of Gunn-type oscillations and the suppression of Bloch oscillations 
is not yet well understood despite theoretical and experimental efforts [3J. 

Although mathematical models at the level of semiclassical kinetic theory go 
back to the 1970s [19], their analysis has been based on simplified reduced 
ordinary differential equations [L41T5] which typically ignore space charge ef- 
fects. Electron transport in a single miniband SL can be described by a kinetic 
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equation coupled to a Poisson equation approximately describing the electric 
potential due to the other electrons [3]. A simple kinetic equation [14J con- 
tains an energy-dissipating collision term of Bhatnagar-Gross-Krook (BGK) 
type PQ and a simple energy-conserving (but momentum-dissipating) collision 
term. This model does not include coupling to the Poisson equation. An im- 
portant point is that the dispersion relation between miniband energy and 
momentum is periodic because this periodicity gives rise to a relation between 
electron drift velocity and electric field which has a maximum value [1] . Then 
the drift velocity decreases as the field increases for large field values (negative 
differential mobility) and this in turn causes the Gunn-type self-sustained cur- 
rent oscillations (SSCO) for appropriate bias and contact boundary conditions 
[I]. These features are absent in the more usual Boltzmann- Poisson systems 
with parabolic band dispersion relations. 

Recently, Bonilla et al [3j have derived a nonlinear drift-diffusion equation 
from the KSS-BGK kinetic model coupled to the Poisson equation, which we 
will call the BGK-Poisson system. They use a Chapman- Enskog perturbation 
method in a particular limit in which the collision terms are of the same or- 
der as the term containing the electric field and dominate all other terms in 
the kinetic equation. Then stable SSCO are obtained by numerically solving 
the drift-diffusion equation with appropriate boundary and initial conditions 
[3]. However, no one has solved numerically the kinetic equation directly and 
shown that self-oscillations are among its solutions or studied the relation be- 
tween these solutions and those of the limiting drift- diffusion equation. These 
are the problems tackled in the present paper and solving them could be a 
step in more precise studies of stable current oscillations in superlattices and 
other low dimensional solid state systems. 

We solve the BGK-Poisson kinetic equation model by means of a deterministic 
weighted particle method that has been used in the past to solve Boltzmann 
equations with non-periodic energy band dispersion relations (20]. Particle 
methods (see a recent one in [10J) are appropriate to study our system of 
equations because their solutions may present large gradients: the electric field 
pulses obtained by simulating the approximate drift- diffusion equations have 
a smooth leading front but a steep trailing back front [I]. The present work 
paves the way to numerically solving interesting problems in nanoelectronics 
and spintronics that are described by related quantum kinetic equations with 
more than one miniband [2]. 



2 The Model 

Our model for electron transport in a single miniband SL is a Boltzmann- 
Poisson system with BGK collision term |1] plus appropriate boundary and 
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initial conditions. The governing equations are: 

dj + v{k)d x f + ^-d k f = -v en (f - f FD ) - ^[/ - /(*, —k, t)}, (1) 

ed 2 x V = ~(n-N D ), F = d x V, (2) 

, I 1 1 



/ FD (fc;/i) = ^%^ln 

7TAT 



(4) 



with x G [0, L] and / periodic in k with period 2tt/1. Here I, L — Nl, N, e, 
/, n, iVo, fee, T, V, — F, m* and — e < are the SL period, the SL length, 
the number of SL periods, the dielectric constant, the one-particle distribu- 
tion function, the 2D electron density, the 2D doping density, the Boltzmann 
constant, the lattice temperature, the electric potential, the electric field, the 
effective mass of the electron, and the electron charge, respectively. We shall 
describe boundary and initial conditions later. 

The first term in the right hand side of Eq. ([T]) represents energy relaxation 
towards a ID effective Fermi-Dirac distribution f FD (k; /x(n)) [3] (local equi- 
librium) due to e.g. phonon scattering. v en is the collision frequency, taken as 
constant for simplicity. Here, /i(n) is the chemical potential that is a function 
of n resulting from solving equation ([3]) when (j3J) is substituted in it. A similar 
BGK model with a Boltzmann local distribution function was proposed by Ig- 
natov and Shashkin [T4"fT5] . The second term in the right hand side of Eq. ([T|) 
accounts for impurity elastic collisions with the constant collision frequency 
Ui mp , which conserve energy but dissipate momentum [T9f3|4"] . Transfer of lat- 
eral momentum due to impurity scattering [T2] is ignored in this model. We 
assume the simple tight-binding miniband dispersion relation, 

e ( k ) = ~2 i 1 ~ cos kl ) =^ v ( k ) = \ ^1 = i^r sin kl > ( 5 ) 

where A is the miniband width. The exact and Fermi-Dirac distribution func- 
tions, / and f FD , have the same electron density n, according to ([3]). The 
latter equation is solved for the chemical potential \x in terms of n, which 
yields the function n{n). When ([1]) is integrated over k, we obtain the charge 
continuity equation, 



d t n + -d x J n = 0, with (6) 

e 
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w/l 

Jn = ^- V(k)f(x,k,t)dk, (7) 
Z7T J 

-it/I 

where J n is the electron current density. 
2.1 Voltage bias condition 



Using the Poisson equation (j2J) to eliminate n, we obtain the following form 
of Ampere's law: 



ed t F + J n = J{t), (8) 

where J(t) is the total current density. The total current density can be ob- 
tained from the voltage bias condition: 

i r 

= - J F(x,t)dx = 0, (9) 



where is the voltage between the two contacts at the end of the SL and 

is an average field. For dc voltage bias, is a fixed constant 0. If we 
integrate (jHJ) over x and use (EJ), we obtain 

i r 

J{t) = -JJ n {x,t)dx. (10) 



2.2 Boundary conditions 



The boundary conditions give the distribution function / on the contacts at 
x = and x = L through the distribution function inside the semiconductor. 
For fixed \k\, there are two possible characteristic curves at a point (x,t): one 
for k > and another one for k < 0. With k < the characteristic curve for 
x — > 0+ and t > is given by the initial condition whereas it is given by the 
distribution function at the contact (x — 0) if k > 0. Then, for i = 0we need 
to specify the distribution function at the contact for k > 0, / + , whereas for 
x = L we need to specify the distribution function at the contact for k < 0, 
f~. Instead of inventing a theory for injecting and collecting contacts, we use a 
top-down approach proposed in Ref . [1] : we know that the following boundary 
conditions appropriately describe current self-oscillations in the drift-diffusion 
equation for the electric field, 
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J n (0,t)=aF(0,t), 
n{L,t) = N D , 



(11) 
(12) 



where a > is the constant contact conductivity and the left hand side of 
Eq. (11 II) is the electron current density. We will use boundary conditions for 
/ such that they become ffTTl) and f[T2"j) when we integrate them according to 
the definitions ([3]) and ([7]) of n and J n respectively: 

f+ = - , /W / v(k)f- dk, (13) 



I J v(k)f(o)dk_^_ 



for x = 0, and 

f _ f® 



(Z/2tt) J\ /(°) d£; 



/ 

V 



/ 

27 



, (14) 



for x = L. Note that the integral over k of f[T3"|) times ev(k) /(2n) yields (TTTT) 
and the integral over of (TH|) times l/(2n) yields f[T2"j) . In these equations, 
is the leading order approximation for the distribution function in the 
Chapman- Enskog method [3~P] : 

oo 

fW(k;n,F)= /fexpfaM), (15) 



where 



AO) = l-ljy/Te rFD 
1 + j V j 



f? D = l -jf D cos{jkl)dk 



F hJv en {v en + ^ imp ) v 

¥ = -ET- F M = — : T e = J 1 + V -. 

F M el V ^en 



Eq. (fl5l) is the solution of ([I]) when we drop the x and t derivatives of / (see 
13J). 

If we use the electric potential V instead of the field F = dV/dx (recall that 
the true electric field is —F), the following boundary conditions for V are 
compatible with (jUJ): 
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L 

V{0, t) = 0, V(L, t)=(t>L = J F{x, t) dx. 

o 



(16) 



2.3 Initial condition 



We select ( TE)|) as our initial condition for the distribution function. The initial 
electric field is assumed to be constant, F(x, 0) = 0, where <fi is the average 
field. If we start from other initial conditions, the evolution of the current and 
other magnitudes are similar to those presented here after about 0.3 ps. 

Recapitulating, the equations governing our model are ([1]) - (0J for the un- 
knowns / and V with initial condition (fT5l) and boundary conditions (|T3|) . 
( TT4l) and ( 1T61) . If we use the field F instead of the electric potential V, the 
voltage bias condition (EJ) for F replaces (TIE]) . 



3 Nondimensional equations 



We use the scales defined in Table 1 to nondimensionalize the Boltzmann- 
BGK-Poisson kinetic equations. These scales are based on the hyperbolic scal- 
ing explained in Ref. [3] 
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Table 1 : Hyperbolic scaling. 
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v m 



Ahr e % (M) ' 

Qfj- (M) = — J cos(jk a ) In [l + exp (m - 5 + 5 cos(k a ] 



dk a , 



where M verifies 



a 



1 = — / In 

2vr 



1 + exp (M - 5 + 5 cos(fc a ; 



dk° 
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with 



a 



m*k B T 
nh 2 N D ' 



Numerical values for these parameters will be given in Section [5j Equations 
(CQ) - have the following nondimensional form 



Al r 

d t «r + ^— sm(k a )d x «r + -F a d k «r 

2hv M V 



f FDa (k a ; n a {n a )) - (l + ^) f a + ^f a (x a , ~k a , t a ) 



dl a V a = d T aF a = n a -l 



(17) 
(18) 



7T 7T 

n a = — [ f a (x a } k a X)dk a = — { f FDa (k a ;fi a (n a ))dk a 
2ix J 2ix J 



(19) 



f FDa {k a - fi a ) = a In [1 + exp (/i a - 5 + 5 cos(fc a ))] 



A 



2k B T 



(20) 



The dimensionless boundary conditions are, for x a = 0: 



/ a + = (3F a - 



f 



a(0) 



Jo" sin {k a ) f a (°) dk 



- J sin(r) f a ~dk a 



(21) 



with 



2-nhaF, 



M 



eAN 



D 



and for x a = L/xq: 



r 



f 



a(0) 



J / t a + AU a 

(l/(2TT))J° 7T f a Wdk a \ 2tt 



1 w 

— [ f a+ dk° 

2vr J J 



(22) 



The boundary conditions for the electric potential V a are 



V a (0, t a ) = 0, V a (L a , t a ) = <p a L a 



F M xq 



(23) 
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The dimensionless initial condition is 



f (°)(^; n «) = £ exp 1 f^ ff Da (n«) (24) 



ff va (n a ) = - / r Wa (A; a ;/i a (n a )) cos(jk a ) dk a 



with x a e [0, L a = L/xq] and / a periodic in fc a with period 2ti. 

Besides the electron current density, J n , it is convenient to calculate the aver- 
age energy E (and its nondimensional version, E a ), defined as E a = E/(kBT): 

E a = S%e{k)f{x,k,t)dk = 5 J^(l - cosk a )f a (x a ,k a ,t a )dk a 
" k B T f% l fl f(x,k,t)dk f^f a (x*,k a ,t«)dk a 



From now on we drop the superscript a. 



4 The Deterministic Weighted Particle Method 



The most widely used numerical method used for solving Boltzmann equations 
is the Monte-Carlo Method [TT] . This stochastic method yields data with a 
lot of numerical noise. The deterministic Weighted Particle Method (WPM) 
is an interesting alternative because it yields the distribution function (and 
therefore its moments: electron density, average energy and current density) 
at each time during the transient regimes with much less noise than the Monte 
Carlo simulation; cf. [2lH[71lo] (a numerical analysis of WPM can be found in 
[S] and in [TB] for the special case of the BGK equation of gas dynamics). 

The WPM relies on a particle description of the distribution function, which 
means that f(x, k,t) is written as a sum of delta functions 

TV 

f(x, k, t) pn^2(Jifi{t)S(x - Xi(t)) ® S{k - fci(t)) 

where u>i, fi(t), Xi(t) and kiit) are, respectively, the (constant) control volume, 
the weight, the position and the wave vector of the ith particle. iV is the 
number of numerical particles. 

In the WPM, the motion of particles is governed by collisionless dynamics, 
whereas the collisions are accounted for by the variation of weights. Large 
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gradients in the solution profile arise from appropriate particles acquiring large 
weights, not by accumulating many particles in the large gradient regions. The 
evolution of the particles is determined by their positions and wave vectors 
which are the characteristic curves of the convective part of the equation. 
Their equations are: 

±kt(t) = ^F,(t), ± Xi (t) = sin (kit)) (26) 

at i] at Invu 



where Fi(t) = F(xi(t),t) denotes the electric field at the instantaneous position 
of the i-th particle. 



The evolution of the weight is given by the ordinary differential equation: 



-£/.<«) = - 

at r] 



1 + 



'invp 



2v P 



v, 



Mt) + f^f (x t (t),-Ht),t) + fnt) 



(27) 



with ff D (t) the Fermi-Dirac distribution evaluated for the i-th particle. 



The system of ordinary differential equations 
a modified (semi-implicit) Euler method: 



(12"7|) is now solved by using 



f? = fr 1 + dt- 

V 



1 + 



v, 



imp 



1v P 



FD,n-l 



2^ e „ 



(28) 



with/r 1 = /(a:r 1 ,-*r 1 ,* B - 1 ) 



e Tpn-1 



K = k?- 1 + dt-Fl . 

V 



x " = % T X + dt 



Al 

2hv 



sin (k\ 



(29) 
(30) 



M 



For stability reasons, we use kf to update x™. The standard Euler method 
would use /c"" 1 to update xf but this would require using unpractically small 
time steps to have a stable scheme. The same problem appears when we employ 
explicit Runge-Kutta or multi-step methods. To select the initial positions 
and wave vectors in the modified Euler method, we build a grid in the domain 
[0, L] x [—Ti, 7r] and choose the values (x°, as the cell centers. The weights 
ff are then chosen according to (I2"lj) . 

The boundary conditions are taken into account as follows: 

• If k™ > 7T, we set kf = kf - 2vr. If kf < -vr, we set kf = kf + 2tt. 
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• If xf > L, we set x™ = x™ — L and = /+. If xf < 0, we set xf = xf + L 
and = f~. Here /+ and /~ are calculated by discretization of the 

integrals in fl2TT) and fl22|) using the composite Simpson's rule on an equally 
spaced mesh K m i with step A/c. 

To calculate Xi, ki and fi at the next time step t n+1 , we need to update 
the electric field and the Fermi-Dirac distribution in the equations for the 
particles. According to Eqs. ([2]) and ([3]), this updating requires an interpolation 
procedure to generate an approximation of the distribution function on a 
regular mesh X m , K m < which is then used to approximate the electric field and 
the chemical potential. To approximate the values of the distribution function 
over the mesh, f mm i, we use the following weighted mean of its values for the 
particles, 

N 

y f n w i , 

t—i J % " ' m,m' 

f n — i=l (Q-n 

J m,m' n / 

V w i , 

L-i ' ' m,m' 



where 



[ \X m -xf\} [ \K m ,-B 

W m ,m' = max 0,1- J ^ *' -max 0,1- ' ^ ' 



and Ax and Ak are the spatial and wave vector steps. 



An approximation for the density ([19]) and average energy f[25|) at the mesh 
points, n(X m ,t n ) « n n m and {k B T) 1 E (X m , t n ) « (ksT) 1 E^, are obtained 
using the composite Simpson's rule and the interpolated values of the distri- 
bution function on the mesh. 

We calculate the nondimensional chemical potential /i by using a Newton- 
Raphson iterative scheme to solve equations (JTJ5J) and (1201 : 

9 (/V-i) ( „ 9) 
9 \Vp-i) 



with 



7T 

g (/I) = n / In [1 + exp (ji — 5 + 5 cos (k))\ dk 

2,71 J 

— 7T 

, . . a f exp (u — 5 + 5 cos (7c)) 

9 ^ ) = -^It 



+ exp (/i — 5 + 5 cos (A;)) 
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The initial guess for \i is obtained by plotting g{p) and selecting a value near its 
zero. g[\x) and g'{p) are evaluated using the composite Simpson's rule. Once we 
know the chemical potential /x, Eq. (j20p provides the Fermi-Dirac distribution 
function at mesh points, f FD {K m i\ /j(n^J), which is then interpolated to get 
the Fermi-Dirac weight function for the particles, ff D,n - 



i/iK,, ( X m+ i — xf \ I K m i + i — k i i Fj) 
Ax J I Ak 



+ 
+ 
+ 



x™ — X m \ ( K m / + i — k™ 



Ax 



Ak 



f FD (K m ,;v(rC +1 )] 

f FD (iW;MO) 



Ax J \ Ak 

\ / ' h n — K 
mil" 1 ! - 1 v m' \ rFD I in 



Ax ) \ Ak 
provided the particle i is in [X m , X m+ i\ x [K m i,K m r+i\. 



+in ' 



(33) 



To compute the electric field at time t n , we use finite differences to discretize 
the Poisson equation on the grid X m : 



V. 



m+l 



2V n + V 



rn—1 



(Ax) 2 

yn _ yn 
" m+l ' m— 1 

2 Ax 



m 



(34) 
(35) 



Here V (0, t n ) = and V (L, t n ) = <\>L as indicated by ((23}. ^ and F™ denote 
our approximations of (X m , t n ) and F (X m , t n ) on the equally spaced mesh 
X m . Finally, the electric field is interpolated at the location of the particle % 



I vi i X-m+x Xj \ fx i X n 

' ^ ( Ax ) /; " + {—AT 



± m+l' 



(36) 



provided the particle i is in [X m ,X m+ i\. 

The total current density J is given by Eq. (flQl) . whose nondimensional version 
is 



±j TV 

J(t) = LJ J sin(k)f(x, k, t) dk 



dx, 



(37) 



in which 

lA 

Anhv M ' 
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We use the composite Simpson rule to approximate J{t n ). 
Summarizing, at each time step t n : 

(1) Calculate the boundary conditions (I2T1) and (1221) with data at time t n ~ l . 

(2) Compute k™ and according to (1281 . (1291) and (!30|) . respectively, by 
using their values at i n_1 . 

(3) Evaluate the distribution function /™ m , at the mesh points (X m ,K m r) 
by the weighted mean fl3TT) . 

(4) Compute the electron density ffl~9]) and nondimensional average energy 
( 125]) at the mesh points. 

(5) Calculate the chemical potential (I3"2"j) and compute the Fermi-Dirac dis- 
tribution f[2"Uj) at the mesh points. 

(6) Interpolate the Fermi-Dirac distribution (|20|) at the mesh points to obtain 
the Fermi-Dirac weight function ff D,n according to (|33|) . 

(7) Compute the electric field at the mesh points by solving the finite differ- 
ence discretization of the Poisson equation, (1341) and (1351) and interpolate 
at the particles according to (13"6T) . 

(8) Calculate the current by evaluation of (13TI) using the composite Simpson's 
rule. 

We have observed that the costlier processes are 2 (computation of /" using 
fl28l) ) and 6 (computation of the Fermi-Dirac weight function ff D ' n ): these 
two processes take about 50% of the overall computation time and they are 
equally costly. After these processes, 3, 5 and 7 have the largest computational 
cost (each takes between 10% and 19% of the overall computation time). 



5 Numerical results 

We have used the parameter values of [9]. Numerical solutions of the non- 
linear drift-diffusion equation derived from the Boltzmann-BGK model show 
that there is a stable stationary state for voltage bias below a certain thresh- 
old. Above this critical voltage, stable self-sustained oscillations of the current 
appear. These oscillations are due to the periodic generation of electric field 
pulses at the injecting contact and their motion towards the receiving con- 
tact. We have observed the same phenomena in our numerical solutions of the 
Boltzmann-BGK kinetic equations. Firstly, we present a typical case of self- 
sustained current oscillations accompanied by the motion and recycling of an 
electric field dipole wave, corresponding to a 157-period 3.64 nm GaAs/0.93 
nm AlAs SL at 14 K, with A = 72 meV, N D = 4.57 x 10 10 cm" 2 , v imp = 2u en = 
18 x 10 12 Hz and dimensionless dc average field = 1 [9]. The constant conduc- 
tivity is 2.5 Q cm -1 and the effective mass is m* = (0. 067 dw + 0.15ds)mo/i, 
where m = 9.109534 x 10~ 31 kg is the electron rest mass. Using these numer- 
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Fig. 1. Total current density versus time plot exhibiting self-sustained oscillations. 
Units are written in parentheses. The oscillation period is 24 ps and the ratio be- 
tween the maximum and the minimum current is 2.6. At the times marked (a) - (f) 
within one oscillation period (30, 36, 42, 48, 54 and 60 ps, respectively), we shall 
depict the electric field profile, the electron density profile, the distribution function 
and the density plots thereof in Figures [21 0E] and [7J respectively. 

ical values, the scales of space, time, velocity, electric field and dimensionless 
chemical potential denned in Section [3] take on the following values: 



x = 15.92 nm, t = 0.23 ps, v M = 68.33 km/s, F M = 22.45 kV/cm, M = 7.11. 



For these parameter values, we consider 140800 particles and a mesh of 440 
grid points for x and 80 points for k. The time step (dt) is 0.002 ps. Figure 
[1] shows the self-oscillations of the current, and Figure [2] the corresponding 
electric field pulse at different times. We observe how the electric field pulses 
are periodically created at the injecting contact x = 0, move to the end of 
the SL and disappear at the receiving contact. In Fig. [21 we have depicted 
the field profiles at the times marked (a) - (e) in FigJU We observe that the 
total current density reaches its maximum value when the electric field pulse is 
about to disappear at the collector. The electric field as a function of time and 
position is shown in Figure [31 both during one oscillation period in Fig. [3]^a) 
and during several periods in Fig. [3](b). The ratio from the maximum to the 
minimum current in FigfJJis 2.6 whereas the same ratio calculated by solving 
the drift-diffusion equation derived in [9] is 2.1 (cf. dashed line in Fig. 1(a) of 
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Fig. 2. Electric field versus position at different times within one period of the 
oscillation. Far from the contacts, at time (c), the electric field pulse is 320 nm wide 
and 139 kV/cm tall. Thus it occupies about 45% of the SL extension. At time (e), 
the electric field has a maximum value of 312 kV/cm. 

[9]). Measured in units of to (which has a different numerical value in [9]), the 
oscillation period is 104.3 in Fig. [1] whereas the drift-diffusion equation yields 
113.8. Comparing Fig. 1(b) of [9j with our Fig. [2j we find that at the time 
(c) the solution of the BGK-Poisson equation produces a pulse far from the 
contacts which is 11 £o wide and 7 Fm tall whereas the drift-diffusion equation 
yields a similar pulse which is 10.7 xq wide and 6.8 tall (cf. dashed line 
in Fig. 1(b) of [9]). Thus the agreement between the simulation of the BGK- 
Poisson system and that of the drift-diffusion equation is very good considering 
the approximations made in the derivation of the latter from the former. 

Figure H] shows the dimensionless electron density. We see the profile during 
several times belonging to one oscillation period as a function of position in 
Fig. 11(a) and as a function of the time and position in Fig. 11(b) . The electron 
density profile corresponding to an electric field pulse is that of a traveling 
dipole wave such that n > 1 behind the peak of the electric field and < n < 1 
ahead of the peak. Comparison with Fig. [3(a) shows that the local maximum 
of the electron density is reached somewhat later than the peak of the electric 
field pulse. 

Figure[5](a) depicts the nondimensional average energy, E /(hsT) as a function 
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Fig. 4. Electron density profiles during one oscillation period. 
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Fig. 5. Nondimensional average energy profiles, E/^bT), at different times of one 
oscillation period. 



of distance at different instants of one oscillation period. The average energy 
profile is pulse-like. Its local maximum is always quite close to the peak of the 
electric field during each oscillation period. Fig. G2h) shows the average energy 
profile as a function of position and time during one oscillation period. 

In Figure El we have depicted snapshots of the distribution function f(x, k, t) 
for different times as marked in Fig. [1] (30 ps, 36 ps, 42 ps, 48 ps, 54 ps, 60 
ps) during one period of the self-oscillations. The structure of the distribution 
function is shown more clearly in the density plots depicted in Fig. [7| for the 
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Fig. 6. (a) - (f ) Distribution function versus position and wave vector at the different 
times of one oscillation period as marked in Fig. (TJ 

same times. The electron density profiles at the these times are shown in 
Figure [SJ We observe that the distribution function has a local maximum at 
location of the peak of electron density. Similarly, / and n have local minima 
at the same positions. The distribution function has a local maximum at a 
positive k (cf. Fig.[7j), and this situation persists from the initial time onwards; 
cf. Figure 

Figure [TUT a) shows the time evolution of the position for six particles, whereas 
Figure [TWb) shows the wave vector vs position for the same particles. The mo- 
tion of the particles is a superposition of an uniform motion and an oscillation 
about it. Comparing Figure ITUlfa) with Fig. [3]^a), we observe that the particle 
positions oscillate with very small amplitudes when the electric field has a 
local maximum at their locations and these amplitudes become larger once 
the pulse has surpassed the particles. In contrast with these great changes in 
oscillation amplitude of the particle positions, the wave vectors of the par- 
ticles oscillate with almost constant amplitudes, as shown in Fig. [TUlfb). [TT1 
and [121 Since the evolution of the particle wave vector is more regular than 
the evolution of the particle position, we can save mesh points on the wave 
vectors. Recalling that the wave vector is a periodic variable, its boundary 
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Fig. 7. (a) - (f) Density plots of the distribution function versus position and wave 
vector at the different times of one oscillation period as marked in Fig. [TJ 




condition is as follows: when one particle goes out of the domain at k = ir, it 
is reintroduced at k = — 7r. This condition can be readily observed in Figures 
dHanddl 
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Fig. 9. Distribution function versus wave vector at t = ps. 




Fig. 10. Time evolution of the position for six particles whose initial dimensionless 
wave vector is —3.126: (a) position vs time, (b) wave vector vs position for the same 
particles. 



6 Convergence of the method 

We have checked the convergence of the method in terms of the number of 
particles N, mesh points M x and and time step dt. Since the Fermi- 
Dirac weights in (1281) and the electric field in (1291) have to be calculated by 
interpolation over mesh points, all these parameters are important for the 
convergence of the calculations over particles and of the calculations over the 
mesh. The physical parameters are the same as in the previous section. 
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Fig. 11. Time evolution of the particle with initial dimensionless position and wave 
vector (13.986, —3.126) during several oscillation periods. 




Fig. 12. Time evolution of the particle with initial dimensionless position and wave 
vector (13.986, —3.126) during one oscillation period. 

Figure [13] shows the evolution of the current density for different number of 
particles N. We have kept the time step and the wave number mesh points 
fixed at the values M k = 80 and dt = 0.008 ps. We observe that we need 
different N for convergence of the calculation depending on the value of M x . 
In Fig. [TBTa). for M x = 440 position mesh points, we have chosen N so that 
N/(M x M k ) takes on the values 1.5, 1.84, 2.25 and 3, whereas in Fig. [T37b). 
M x = 360 and N is chosen so that N/ (M x M k ) takes on the values 1.5, 2.25 and 
3. We observe that for dt = 0.008 ps and N/(M x M k ) > 2.25 the results do not 
change if we increase the number of particles. In particular, for N = 64800, 
N/(M x M k ) = 1.84 if M x = 440 and N/(M x M k ) = 2.25 if M x = 360. In 
the first case shown in Fig. [TBTa). we need more particles for the method to 
converge whereas in the second case, Fig. [TBTb) shows that we do not improve 
our results by increasing N. The convergence range of N/ (M x M k ) depends 
slightly on the time step dt: if dt = 0.002 ps, M x = 440, M k = 80, our 
calculations yield indistinguishable curves J(t) for Nj \M x M k ) = 3, 4, but 
not for N/(M x M k ) = 2.25. Thus we have found that it is advisable to select 
iV so that 3 < N/(M x M k ) < 4.5: numerical results are indistinguishable 
when N/ (M x M k ) > 3 and the computational cost is not very large if we keep 
N/(M x M k ) < 4.5. 

Figures PT41 and IT51 show the evolution of the total current density in simulations 
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13. Current versus time for different number of particles N when = 80 
0.008 ps. (a) M x = 440 and (b) M x = 360. 
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Fig. 14. Current versus time for different number of wave vector mesh points when 
M x = 520, N = 200000 and dt = 0.008 ps. 
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Fig. 15. Current versus time for different number of position mesh points when 
M fc = 80, N = 480000 and dt = 0.008 ps. 
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Fig. 16. Current versus time for different time steps when N = 90000, M x = 260 
and M k = 80. 

having a different number of mesh points when the time step is dt = 0.008 
ps. In Fig. [HI different wave vector mesh points are considered for M x = 520 
and N = 200000. We can check that we do not need to have a very fine grid 
in k because we obtain the same results with Mk = 60 and Mk = 80. In Fig. 
[T5l different numbers of position mesh points are considered for Mk = 80 and 
N = 480000. For M x = 920 and larger our numerical curves overlap. 

Lastly, Figure [16] shows the evolution of the total current density for different 
time steps in simulations with 90000 particles and 260 mesh grid points for the 
position and 80 for the wave vector. We observe that our results are similar 
for time steps dt = 8 x 10~ 4 ps and smaller. 

Figures [14] to [16] show that the shape of J(t) is similar for different mesh 
points and time steps: the device behavior is qualitatively correct even if we 
take fewer mesh points or larger time steps than needed to attain a numerically 
precise current vs time graph. Smaller Mk, M x and larger dt result in slightly 
smaller oscillation periods and slightly larger oscillation amplitudes. 

Our numerical simulations have been carried out using a Matlab code in a 
computer with a Genuine Intel(R) CPU T2050 @ 1.60GHz processor with a 
1595 MHz speed. Several computation times for time steps dt of 0.008 and 
0.002 ps and 10000 time steps are shown in Table 2. Clearly, the time the 
computer takes to calculate one time step dt decreases as the number of par- 
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tides, M x or Mk decrease. Except for the last row in Table 2, all rows satisfy 
N/ (M x Mk) > 2.25, and the corresponding particle numbers and x and k mesh 
points produce accurate enough results. 



N 


M x 


M k 


AT 

M x M k 


dt (ps) 


(C.T.)/# steps (seconds) 


480000 


920 


80 


6.52 


0.008 


2.64 


480000 


360 


80 


16.67 


0.008 


2.12 


200000 


520 


80 


4.81 


0.008 


1.19 


200000 


520 


32 


12.02 


0.008 


1.15 


140800 


440 


80 


4.00 


0.002 


0.87 


105600 


440 


80 


3.00 


0.002 


0.73 


79200 


440 


80 


2.25 


0.002 


0.65 


52800 


440 


80 


1.50 


0.002 


0.53 



Table 2 : Computer time (C.T.). 



7 Conclusion 

We have proposed a deterministic weighted particle method to numerically 
solve for the first time the semiclassical Boltzmann-BGK-Poisson system of 
equations with periodic miniband energy dispersion relation. This system de- 
scribes vertical electron transport in a GaAs/AlAs superlattice under dc volt- 
age bias conditions. When using appropriate values for the injecting contact 
conductivity and voltage, we find a stable self-sustained oscillation of the cur- 
rent through the structure which corresponds to periodic nucleation of electric 
field pulses at the injecting contact that then move to the receiving contact. 
The pulses have a large electron density on their trailing edges which implies 
large gradients of the electric field there. These gradients are well resolved by 
particles having large weights there, which is one of the advantages of using 
the weighted particle numerical method. Our results agree with experimental 
observations [T3|4] and confirm the validity of the Chapman-Enskog pertur- 
bation method used to derive a drift-diffusion equation for high electric fields 
[3]. In fact, the electric field profile and the total current density obtained by 
numerically solving the the drift- diffusion equation [9] agree very well with 
the numerical solution of the kinetic equations obtained in the present work. 
Having solved the kinetic equations directly, we can obtain the evolution of 
the distribution function and its relevant moments such as electron density, 
current density and average energy. The present work paves the way to numer- 
ically solving interesting problems in nanoelectronics and spintronics that are 
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described by related quantum kinetic equations with more than one miniband 
0. 
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